Large-Scale Simulations of Diffusion-Limited n-Species Annihilation 
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We present results from computer simulations for diffusion-limited n-species annihilation, Ai + 
Aj — » (i, j = 1, 2, . . . , n; i j), on the line, for lattices of up to 2 28 sites, and where the process 
proceeds to completion (no further reactions possible), involving up to 10 15 time steps. These 
enormous simulations are made possible by the renormalized reaction-cell method (RRC). Our 
results suggest that the concentration decay exponent for n species is a(n) = (n — l)/2n instead of 
(2n — 3)/(4n — 4), as previously believed, and are in agreement with recent theoretical arguments [IcH|. 
We also propose a scaling relation for A, the correction-to-scaling exponent for the concentration 
decay; c(t) ~ t~ a (A + Bt~ A ). 

PACS numbers: 05.40.-a, 82.20.Wt 
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I. INTRODUCTION 

Diffusion-limited reactions have attracted much inter- 
est in recent years P, Q] . The kinetics of such systems is 
dominated by local fluctuations in the concentration of 
the reactants, thus posing a formidable problem which 
has not yet been solved: there exists no comprehensive 
theoretical approach for the analysis of diffusion-limited 
processes. 

Few select models are amenable to exact analysis. 
These include one-species annihilation, A + A — > 0, and 
two-species annihilation A + B — > (see [3J and refer- 
ences therein). In one dimension, the particle density 
for one-species annihilation decays as c(t) ~ t~ 1 / 2 , while 
for two-species annihilation (with equal initial concentra- 
tions and same diffusion constants for the two species) 
c(t) ~ i" 1 / 4 . In either case the result is strikingly dif- 
ferent from the mean-field kinetics of the corresponding 
reaction-limited process, c(t) ~ 1/t. To bridge the gap 
between these disparate behaviors, ben-Avraham and 
Redner Q proposed the n-species annihilation model, 
where particles belonging to the n species Ai , A%, . . . , A n 
diffuse on the line and react immediately uppon en- 
counter, according to the scheme: 



Ai + Aj 



0. 



,.7 = 1,2, 



(1) 



For n = 2 we recover two-species annihilation, while in 
the limit n — > oo encounters between like-particles are 
improbable and the model is equivalent to one-species 
annihilation. For intermediate values of n, one expects 
c(t) ~ 

In it was proposed, following a heuristic scaling ar- 
gument and treating fluctuations via the Van Kampen 
fi-expansion [f|, that a(n) = (2n — 3) /(An — 4). This was 
supported by numerical simulations of lattices of typi- 
cally 10 6 sites, and up to 10 6 time steps. (In one time 
step, all of the particles in the system move one lattice 



spacing each, on average.) Recently, we have conducted 
extensive numerical simulations , followingthe method 
of Renormalized Reaction-Cells (RRC) U&M The sys- 
tems involved are up to 2 28 w 2.7 x 10 s sites long, and 
the processes were simulated to completion (until no fur- 
ther reactions are possible), for up to 10 15 time steps. 
The new data leads us to the conjecture that a(n) = 
(n — l)/2n 6J. We also find a correction to the main 
decay mode, of the form c(t) ~ t~ a( - n \A + Bt~ A ^), 
A(n) = l/2n. The same results were found, inde- 
pendently (and unbeknownst to us), by Deloubriere et 
al., [ljj. I n their theoretical derivation, they consider 
a simplified version of n-species annihilation, where do- 
mains of alternating species loose particles to reactions at 
one and the same rate, in a synchronous fashion. The ap- 
proximation is more than reasonable, yet it does not rig- 
orously apply to the original model, and analysis of cor- 
rections is certainly beyond its scope. Moreover, the sim- 
ulations in ^3] are comparable in size to those in 0. In 
what follows, we report the results of our large-scale sim- 
ulations, which strongly support the conclusions of jToj | . 
We also propose a scaling relation for the correction expo- 
nent A for n-species annihilation, and possibly for other 
reaction models where particles segregate into distinct 
domains. 



II. SCALING 

As is well known, local fluctuations in the concen- 
trations of the different species drive the kinetics of n- 
species annihilation 0, 0, 3 • An initially random homo- 
geneous distribution of the particles evolves into a contin- 
uously growing mosaic of domains of alternating surviv- 
ing species. Two lengthscales characterize the emerging 
distribution and dominate the system: the inter-domain 
distance — the distance between the last particle in a 
domain and the first particle in the domain next to it — 
^ab(^), and the domain length, £(t) These quantities 
grow with time as 
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t(t)~tP, £AB(t)~f>. (2) 



2 



Once domains form, reactions might take place only 
at the domain boundaries, and particles have to dif- 
fuse across the domain gap £ab to react with other 
species. This takes a typical time of At ~ £ AB /D, 
where D is the diffusion constant. The change in par- 
ticle concentration during time At equals the total num- 
ber of domain boundaries divided by the lattice size L; 
Ac - -{L/£)/L = -1/L Thus, 



On substituting the relations J5J) and c{t) ~ t Q , we de- 
rive the scaling rule 

2<y + j3-a = l, (4) 

Due to the underlying transport mechanism, we expect 
that domains grow diffusively, as £ ~ t 1 ^ 2 , so f3 = 1/2, 
and in effect there is only one independent exponent: 
2a — 7 = 1/2. The general scaling form holds also for 
two-species annihilation in the presence of drift (and with 
ard-core repulsion between like species), where a = 1/3, 
13 = 7/12, and 7 = 3/8 @. 



III. SIMULATION RESULTS 

The n-species annihilation process is simulated as fol- 
lows. The sites of a one-dimensional lattice are cither 
empty or occupied by a particle (of one of the n species) . 
Periodic boundary conditions are imposed, so the lattice 
is effectively a ring. At each Monte Carlo step a particle 
is chosen randomly and is moved to the nearest site to 
its right or left, with equal probabilities. If the target 
site is occupied by a particle of a different species, then 
both particles are removed from the system, mimicking 
the reaction 0J. If the target site is occupied by a parti- 
cle of the same species, then the move is disallowed and 
it does not take place. Regardless of the outcome, time is 
incremented by 1/N(t), where N(t) is the total number 
of extant particles. 

As the simulation proceeds, the particle concentration 
declines and the typical distance between particles in- 
creases. The time spent on simulating the diffusive mo- 
tion of the particles until they interact grows even faster, 
as the square of the distance between them. Because 
of that, computer simulations are limited to relatively 
short times. This problem is overcome by the RRC 
method 0,H H. 

In the RRC method the particles occupy cells, rather 
than sites. Each time the concentration halves, the cells 
are renormalized: every two cells are merged into one, 
and time is renormalized accordingly. The typical time 
required to diffuse out of a renormalized cell twice as 
large as that of the previous generation is four times 
longer. Thus, physical time is simulated faster with each 
renormalization step and the process can be simulated to 



completion. Other details for the implementation of the 
RRC method are discussed in 

Simulations were performed on DEC Alpha processors 
running Linux. Since each lattice site requires 6 bytes 
(for species, number of particles, and a pointer to a list 
of populated sites that is used for fast selection at each 
Monte Carlo step), with 2 Gigabytes (2 31 bytes) memory 
we were able to simulate lattices of up to 2 28 sites. The 
compiler was given special #pragma pack(l) instructions 
to circumvent word alignment (which would allocate 32 
bytes for our 6-byte site). 

To test the technique, we have simulated the cases of 
n = 2 and n = 3 on lattices of 2 16 = 65, 536 sites, in both 
the RRC and the traditional simulation method. These 
lattices are small enough to enable the simulation of the 
process by the traditional method to completion. On the 
other hand, the system is large enough to let us exam- 
ine the effect of the renormalizations: with 2 16 sites and 
c(0) = 1/16 the RRC method requires 12 renormaliza- 
tions. In Fig.^, we compare the particle concentration as 
obtained by the two methods. In Fig.^ 3 we plot the local 
slope of the curves, the exponent a(t). The renormaliza- 
tions are discernable only in this second, more stringent 
test, but the overall agreement is excellent. Similar re- 
sults were obtained for the domain size and the distance 
between domains. 
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FIG. 1: Comparison between the traditional (circles) and 
RRC (solid line) simulation methods. Plotted is the num- 
ber of surviving particles N(t) (top) and the local slope 
a(t) = -d\nN(t)/dlat (bottom) for 2 16 -site lattices. 
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Having gained some confidence in the RRC method, 
we proceed to larger simulations. In Fig. [21 we show the 
surviving number of particles, N(t), at time t, for n = 3, 
4, and 5, and several lattice sizes. In Fig. El we plot the 
local decay exponent a(t) for our largest simulations of 
n = 3. The maximum of the curve at t rs 10 4 agrees 
with the earlier prediction that a = 3/8 4]. (Indeed, 
simulations in |4| yielded a somewhat smaller value than 
the theoretical 3/8, in perfect agreement with current 
results.) However, a(t) is seen to diminish with time, 
suggesting a long-time asymptotic limit of a w 1/3. This 
limiting value is confirmed in the data collapse (especially 
at long times) of Fig.^J where we plot t a c(t) vs. t^/L for 
various system sizes, and a=l/3,/3=l/2. Independent 
measurements show that (3 = 1/2, as assumed, to within 
2%, and the data collapse of Fig. ^deteriorates with other 
choices for the values of a and (3. 
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FIG. 4: Scaling of concentration, c(t) — t~ a p(t' 3 /L), for 3- 
species annihilation. The best data collapse (at late times) is 
obtained for a — 1/3 and /3 = 1/2. 
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FIG. 2: Concentration decay for n = 3, 4, and 5-species an- 
nihilation. Plotted is the number of surviving particles, N(t), 
for system sizes L = 2 16 , 2 20 , 2 24 , and 2 28 (bottom to top). 
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FIG. 3: Local decay exponent for 3-species annihilation. 

We have analyzed in this fashion n = 3, 4, and 5- 
species annihilation, and measured the exponents a, (3, 
and 7. Our results are summarized in Table [I] In all 
cases, the scaling relation Q seems to hold, and (3—1/2 
to within numerical errors. Looking for a simple expres- 
sion of these results, that would have the appropriate 
limits for the known cases of n — 2 (two-species annihi- 



lation) and n — ► 00 (one-species annihilation), we were 
led to the conjecture p 

a = —> P = r 7 = ^T' (5) 

a result derived independently by Deloubriere et al. 01 



n q 


n—l 

2n 


P 


7 


•ln-L 

in 


3 0.33(1) 

4 0.39(2) 

5 0.42(2) 


0.333 
0.375 
0.400 


0.50(1) 
0.50(1) 
0.50(1) 


0.42(2) 
0.44(1) 
0.47(2) 


0.417 
0.434 
0.450 



TABLE I: Exponents a, f3 and 7 



Finally, let us address the issue of corrections to scaling 
of the concentration decay. We look for corrections of the 
form 

c(t) ~t- a {A + Bt- A ) , (6) 

where A and B are constants. Our strategy consists of 
performing a least squares linear fit of A+Bt~ A to t a c(t), 
for different powers A, and searching for the value of 
A which minimizes the error. The scaling form @ is 
expected to work only after the asymptotic regime sets 
in, and before finite-size effects begin, and the sticky part 
of our procedure is deciding which times demarcate this 
region. Experimenting with different choices gives us a 
feel for the errors involved. In Fig. we show best fits 
for the region t = 10 6 - 10 12 , for n = 3 on a L = 2 28 
lattice, where our data is most reliable. The results are 
most compatible with A = 1/6 (for n = 3). Similar 
tests for other values of n lead us to the conjecture that 
A(n) = l/2n. 

The correction exponent can be understood by a 
simple-minded argument. In deriving @ we have as- 
sumed that the typical distance between reacting parti- 
cles, at the edges of adjacent domains, is £ab- While this 
is correct, we note that, had the distribution of particles 
been homogeneous, the distance between reacting pairs 
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would be typically Iaa ~ L/c ~ i Q , quite different from 
the assumed €ab ~ Using At ~ £ AA /D in J3J), instead 
of £\ B /D, yields a faster decay; c ~ t _ ( 1 ~ /3 ). Diffusion 
provides a natural drive toward a homogeneous distri- 
bution, and so it is conceivable that this faster mode of 
decay is manifested as a correction to the main behavior, 
c ~ t~ a . It follows from iJBJ that the correction exponent 
is 

A=l-0-a=± (7) 

where the last equality applies to n-species annihilation, 
provided that the conjecture JSJ holds. The more general 
relation works well for two-species annihilation with drift, 
where a = 1/3, /3 = 7/12, and A = 1/12 @. 
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FIG. 5: Corrections to scaling. Simulation results (circles) 
are best fit by Eq. (JSJ, with A = 1/6 (solid line). 
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IV. SUMMARY AND DISCUSSION 



We have presented large scale simulation results for 
diffusion-limited n-species annihilation, in one dimen- 
sion, using the RRC method. Our simulations contra- 
dict previous work and are in favor of new theoretical 
arguments advanced by Deloubriere et al. [13 • We have 
also provided a new scaling relation for the correction-to- 
scaling exponent A, valid for diffusion-limited reactions 
in one dimension, where the particles segregate into al- 
ternating domains. The corrections to the main decay 
mode are large, and explain the failure of 4] to obtain 
the correct asymptotic behavior with the size of simula- 
tions employed at that time. An important conclusion 
to be drawn is that predicting asymptotic behavior from 
the typical size of simulations used commonly in the field 
is dangerous. More advanced techniques and larger sim- 
ulations seem to be imperative. 
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